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Abstract 

In this paper, we consider the numerical solution of poroelasticity problems that are of Biot type 
and develop a general algorithm for solving coupled systems. We discuss the challenges associated with 
mechanics and flow problems in heterogeneous media. The two primary issues being the multiscale nature 
of the media and the solutions of the fluid and mechanics variables traditionally developed with separate 
grids and methods. For the numerical solution we develop and implement a Generalized Multiscale 
Finite Element Method (GMsFEM) that solves problem on a coarse grid by constructing local multiscale 
basis functions. The procedure begins with construction of multiscale bases for both displacement and 
pressure in each coarse block. Using a snapshot space and local spectral problems, we construct a basis 
of reduced dimension. Finally, after multiplying by a multiscale partitions of unity, the multiscale basis 
is constructed in the offline phase and the coarse grid problem then can be solved for arbitrary forcing 
and boundary conditions. We implement this algorithm on two heterogenous media and compute error 
between the multiscale solution with the fine-scale solutions. Randomized oversampling and forcing 
strategies are also tested. 


1 Introduction 

Problems of mechanics and flow in porous media have wide ranging applications in many areas of science 
and engineering. Particularly in geomechanical modeling and its applications to reservoir engineering for 
enhanced production and environmental safety due to overburden subsidence and compaction mm- ° ne 
of the key challenges is the multiscale nature of the geomechanical problems. Heterogeneity of reservoir 
properties should be accurately accounted in the geomechanical model, and this requires a high resolution 
solve that adds many degrees of freedom that can be computationally costly. Moroever, there are disparate 
scales between the often relatively thin reservoir structure and the large overburden surrounding the reservoir 
that adds more complexity to the simulation. Therefore, we propose a multiscale method to attempt overcome 
some of these challenges. 

The basic mathematical structure of the poroelasticity models are usually coupled equations for pressure 
and displacements known as Biot type models m- For pressure, or flow equations, we have the parabolic 
equation Darcy equation with a time dependent coupling to volumetric strain. The stress equation is the 
quasi-static elasticity equations with a coupling to the pressure gradients as a forcing. Poroelastic models of 
this type have been explored in the petroleum engineering literature in the context of geomechanics for some 
time pi on m m to name just a few. There are noted issues that arise. The first being heterogeneities of 
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the reservoir and surrounding media add many complications to the effective simulation due to complexity of 
scales. Moreover, development of flow and mechanics simulation were often considered separately. Progress 
was made on this problem by considering various coupling strategies m- However, in the instance that the 
physics is not well understood a fully coupled scheme may be desired. This separation of development from 
flow and mechanics methods adds the complication of the computational grids not being the same in each 
regime. Some effort has been made in the improvements of gridding techniques between geomechanical and 
flow calculations [20] and references therein. 

As briefly noted before, typically for numerical solution of such coupled systems time splitting schemes 
are often used. Various splitting techniques for poroelastic equations have been explored and analyzed in 
the context of reservoir geomechanics in mm- Also, in the context of poroelasticty and thermoelastic 
equations, various splitting techniques have been analyzed and implemented [15] . The primary splitting 
techniques are the undrained, fixed-stress, and fully implicit. Due to observed better errors, we will primarily 
consider the less computationally costly fixed-stress splitting and the more robust, yet with a loss in some 
matrix sparsity, fully implicit coupled approach. 

Once the equations have been split in time we wish to resolve in space and will utilize a multiscale method. 
There are many very effective multiscale frameworks that have been developed in recent years. There are 
rigorous approaches based on homogenization of partial differential equations, where effective equations 
are derived based fine-scale equations at the microstructure level ei m- However, these approaches may 
have limited computational use and more practical multiscale methods are used. Examples include the 
Heterogeneous Multiscale Method (HMM), where macro-scale equations on coarse-grids are solved while the 
effective coefficients on the fine-scale are resolved at each coarse grid nodes [ 22 ., [23]. An approach based 
on the Variational Multiscale Method (see [21]), where coarse-grid quasi-interpolation operators are used to 
build an orthogonal splitting into a multiscale space and a fine-scale space [25] • Fine-scale space corrections 
are then localized to create a computationally efficient scheme. In this paper, we will use the Generalized 
Multiscale Finite Element Method framework, which is a generalization of the multiscale finite element 
method [2] . 

To efficiently solve these splitting schemes and overcome some of the challenges of heterogenous reservoir 
properties and gridding issues between mechanics and flow, we will develop a Generalized Multiscale Finite 
Element Method (GMsFEM) pQ. Our GMsFEM has the advantage of being able to capture small scale 
features from the heterogeneities into coarse-grid basis functions and offline spaces, as well as having a 
unified computational grids for both mechanics and flow solves. The offline multiscale basis construction 
may proceed in both fluid and mechanics in parallel and both constructions are comparable. We proceed by 
first generating a coarse-grid and in each grid block a local static problem with varying boundary conditions 
is solved to construct the snapshot spaces. We then perform a dimension reduction of the snapshot space 
by solving auxiliary eigenvalue problems. Taking the corresponding smallest eigenpairs, and multiplying by 
a multiscale partition of unity we are able to construct our offline basis. In this greatly reduced dimension 
offline basis, the online solutions may be calculated for pressure and displacements for any viable boundary 
condition or forcing. 

The work is organized as follows. In Section 2 we provide the mathematical background of the poroe- 
lasticity problem. We will introduce the Biot type model and highlight where the heterogeneities primarily 
occur. In our formulation, the computational domain will be entirely inside of the fluid filled, or reservoir, 
region. However, coupling to regimes of pure elasticity to model the overburden are of course possible. In 
Section 3, to outline the difficulties in full direct numerical simulation we introduce the fine-scale discretiza¬ 
tions using coupled and splitted schemes. Once we split the porooelastic system we will be able to apply 
our multiscale method. In Section 4, we present our GMsFEM algorithm and outline its construction pro¬ 
cedure. We will use the offline multiscale basis functions to calculate accurately pressure and displacements, 
at a reduced dimension and computational cost in the online phase. Finally, numerical implementations 
are presented in Section 5. Using the GMsFEM, we compare the multiscale solution to fine-scale solutions 
and give error estimates. We will present two different examples with varying coefficients. Additionally, we 
will implement and discuss different strategies with oversampling and randomized forcings to construct the 
multiscale spaces. 
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2 Problem formulation 


We denote our computational domain B C to be a bounded Lipschitz region. We consider linear poroe- 
lasticity problem where we wish to find a pressure p and displacements u satisfying 


9 div it 


a 


dt 


— div a(u) + a grad(p) 
1 dp (k \ 

+ s a _ d "'ly 


= 0 in fl, 

= / in fi, 


(la) 

(lb) 


with initial condition for pressure p(x,0) = po. We write the boundary of the domain into four sections 
311 = r\ ur 2 = r 3 ur 4 . We suppose the following boundary conditions on each portion 


an = 0, igTi, u = ui, x £ T 2 , 


and 

— w == o, x g r 3 , p = pi, x e r 4 . 

v on 

Here the primary sources of the heterogeneities in the physical properties arise from a, the stress tensor and 
k, the permeability. We denote M to be the Biot modulus, v is the fluid viscosity, and a is the Biot-Willis 
fluid-solid coupling coefficient. Here, / is a source term representing injection or production processes and 
n is the unit normal to the boundary. Body forces, such as gravity, are neglected. In the case of a linear 
elastic stress-strain constitutive relation we have that the stress tensor and symmetric strain gradient may 
be expressed as 

a(u) = 2 pe(u) + A div(u) X, e{u) = - (grad u + grad u T ) , 


where p, A are Lame coefficients, X is the identity tensor. In the case where the media has heterogeneous 
material properties the coefficients p and A may be highly variable. 

The above poroelasticity problem (la), assuming a linear elastic stress-strain relation, can be written in 
operator matrix form: 


Au + aGp = 0, 


-f- (Sp + aDu ) + Bp 
dt 


/, 


( 2 ) 

(3) 


where 


Av = —pX/ 2 v — (A + p) grad div v, Bp = — div ( — gradp 


and G and D are gradient and divergence operators and S = jjX. 


3 Fine-Scale Discretization 

We will now present splitting methods for the above system in the context of solving the fine-scale approx¬ 
imation. This will highlight the areas where we would like to utilize a multiscale method when solving in 
the spatial variables due to the degrees of freedom required in resolving the system. For approximating the 
numerical solution to ([Tj) on fine-scale grid we use a standard finite element method. We begin by giving the 
corresponding variational form of the continuous problem written as 

a(u, v ) + g(p , v) = 0, for all v £ V, (4) 

d +C +b ( p ' q ' > for all q £ Q. (5) 
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for u £ V, p £ Q where 


V = {v £ [i? 1 (fl)] d : v(x) = ui,x £ r 2 }, Q = {q £ ff 1 (Q) : q(x) = Pi,x £ r 4 }, 
and the test spaces with homogeneous boundary conditions are given by 

V = {v £ [ff 1 (f2)] d : v{x) = 0, x £ r 2 }, Q = {q £ H 1 (f2) : q{x) = 0,cc £ r 4 }. 


Here for bilinear and linear forms we have define 


a ( u , v ) = (r(u) v dx, b(p,q) = J^^gradp,gradqj dx, c{p,q) = £ ^ pqdx, 


and 


g{p,v) = / ck(gradp, v)dx, d{u,q) = / adivuqdx, ( f,q)= / fqdx. 

J r2 »fi o 

Here (•, ■) under the integrand denotes the standard inner product. In Section [5j we will discretize the spaces 
using a fine-scale standard FEM and denote them Vh,Qh and Vh,Qh, h being the fine-grid size. The FEM 
using these spaces will serve as a reference solution for our GMsFEM outlined in Section [3] 

To solve the above system we first discretize in time. This discretization leads to several possible couplings 
between time-steps and the two equations of prorelasticity. We proceed by giving the coupled and so-called 
fixed-stress splitting mm- The standard fully implicit finite-difference scheme, or coupled scheme, can be 
used for the tinre-discretization and is given by 


d 


^u n+1 -u n 



a(u n+1 ,v)+g(p n+1 ,v) = 0, 
+ C ( y ~ -^“>9) +b(p n+1 ,q) = (/, q), 


(6a) 

(6b) 


with u n = u(x,t n ), p n = p(x,t n ), where t n = nr, n = 0,1,...,My, Mtt = T and r > 0. For time 
discretization we can apply many different splitting techniques which often occur in the literature. 

Another we shall consider here is the fixed-stress splitting scheme 


d 





a(u n+1 ,v)+g(p n+ \v) = 0, 
+ s ^ +b(p n+1 ,q) = l(q), 


(7a) 

(7b) 


where the variational form is re-written with 


s(p,q) 




pqdx, 




a 2 p n — p n 1 \ 
d^-dr T ) 


q dx 


and Kdr is the drained modulus 


d^dr 


E{ 1 - u p ) 

(i - 2v p )(i + v P y 


where u p is the Poisson ratio and E is the elastic modulus. When we utlize the fixed-stress splitting scheme, 
first we solve pressure equation for p n+1 given data at the previous time-steps. Then, passing this new 
pressure information, we return to the quasi-static stress equation and calculate displacements at u n+1 . 


4 GMsFEM for Poroelasticity 

In the GMsFEM presented here, we will focus on the development in the fixed-stress splitting Q. We will 
however give numerical examples from both coupling strategies. The fixed-stress splitting decouples the 
flow and mechanics equations. We will first present the offline multiscale basis construction in the fluid or 
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pressure solve then its construction in the mechanics or displacement calculation step. In this algorithm, due 
to the heterogeneities arising primarily from the permeability k and the stress tensor <r(u), we will solve local 
problems in each of the relevant portions of the variational form to construct the offline multiscale spaces. 
We now outline the general procedure of the GMsFEM algorithm. 

The overall fine-scale model equations will be solved on a fine-grid using spaces 14, Qh and Vh,Qh, and 
will act as our reference solutions. Once the fine-grid is established we must introduce the concepts of coarse- 
grids and their relationships. To this end, let T H be a standard conforming partition of the computational 
domain fl into finite elements. We refer to this partition as the coarse-grid and assume that each coarse 
element is partitioned into a connected union of fine grid blocks. The fine grid partition will be denoted by 
T , and is by definition a refinement of the coarse grid T H ■ We use , where N is the number of 

coarse nodes, to denote the vertices of the coarse mesh T H , and define the neighborhood of the node Xi by 

^i = \J{ K 3^T a \ XiCKj}. 

3 

See Figure[l]for an illustration of neighborhoods and elements subordinated to the coarse discretization. We 


T h (Coarse Grid) 



Figure 1: Illustration of a coarse neighborhood and coarse element 

emphasize that the use of Wj is to denote a coarse neighborhood, and we use I\ to denote a coarse element 
throughout the paper. 

Boadly speaking, the GMsFEM algorithm consist of several steps: 

• Step 1: Generate the coarse-grid, T H . 

• Step 2: Construct the snapshot space, used to compute an offline space, by solving many local 
problems on the fine-grid. 

• Step 3 : Construct a small dimensional offline space by performing dimension reduction in the space 
of local snapshots. 

• Step 4: Use small dimensional offline space to find the solution of a coarse-grid problem for any force 
term and/or boundary condition. 
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As noted previously, because coupled system of equations for poroelasticity problems can be solved using 
splitting scheme, we can construct multisclate basis functions for pressure and displacements separately. We 
begin by considering the pressure solve, then, the displacement solve. 


4.1 Pressure Solve 


Recall, for the numerical solution of pressure equation on coarse grid we consider the continuous Galerkin 
(CG) formulation (7b I given by 

n n+1 


p 


■p 


q) + b(p n +\q) = l(q)~d 


for all q £ QoS, 


( 8 ) 


where Qoff is used to denote the space spanned by multiscale basis functions > each of which is supported 
in lo 1 . The index k represents the numbering of these multiscale basis functions. We will now show how to 
construct the offline multiscale space Q 0 g. In turn, the CG solution of the form 


KM) = '52p l k{t)i’k i ( x )> 

i,k 


will be sought. 

We begin by construction of a snapshot space E s “ ap . We use harmonic extensions 

b(ipj ,ana ‘ p ,q) =0 in w, 

Vf-nap = 5*(z) on du>. 

Here Sj(x) are defined by 6j(x) = Sj Vj, k £ where Jh{^i) denotes the fine-grid boundary node on 

du>i. For simplicity, we will omit the index i when there is no ambiguity. 

Let l, be the number of functions in the snapshot space in the region w, and define 

snap = span{^ nap : 1 <j<h}, 

for each coarse subdomain w. We denote the corresponding matrix of snapshot functions to be 

^ap=[VC^---,C aP ]- 

To construct the offline space Q 0 g, we perform a dimension reduction of the space of snapshots by using an 
auxiliary spectral decomposition. More precisely, we solve the eigenvalue problem in the space of snapshots: 

B oS ^f = A fM oS ^f, (10) 


where 


''snap/ ''snap 

where B and M denote fine scale matrices 


B° S = (RLapf^nap, M° S = 


''snap 7 


grad (j>i , grad <pj ) dx, M,^ = / —ifii^jdx. 


Bij — 


Here, <j>i are fine-scale basis functions. 

We then choose the smallest N“g P eigenvalues from Eq. (10) and form the corresponding eigenvectors in 
the space of snapshots by setting 

r k s = J2^ s j 

i=i 
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for k = 1.... ,7V"/, where are the coordinates of the vector We denote the span of this reduced 
space as 

For construction of the offline space, to ensure the functions we construct form an H 1 conforming basis, 
we define multiscale partition of unity functions Xi 


b{Xi,Q) = 0 in K, 
Xi = 9i on c)K. 


( 11 ) 


for all K G w. Here g t is a continuous on K and is linear on each edge of dK. We could choose g t to also 
be selected shape function, Neumann conditions, or boundary conditions on larger domains in the context 
of oversampling. 

Finally, we multiply the partition of unity functions by the eigenfunctions in the offline space to 
construct the resulting basis functions 

ipi,k = l° r 1 < * < N and 1 < k < N^’ p , (12) 


where N^g P denotes the number of offline eigenvectors that are chosen for each coarse node i. We note that 
the construction in Eq. (12) yields continuous basis functions due to the multiplication of offline eigenvectors 
with the initial (continuous) partition of unity. Next, we define the continuous Galerkin spectral multiscale 
space as 


QoS = span {ipi t k '■ 1 <i < N and 1 < k < N^’ p }. 


(13) 


Using a single index notation, we may write Q 0 s = spanwhere Nj? = IVos’ P denotes the total 

number of basis functions in the spaces for i = 1 ,,N. 

Denote the matrix 

RP=ty V>iv?] T > 

where ^ are used to denote the nodal values of each basis function defined on the fine grid. Then, the 
variational form in (18]) yields the following linear algebraic system 


QcPc +1 = y cP : , 


(14) 


where 

Qc = RP ^ + Fc = RPFp - 

Here, F p being the operator corresponding right hand side data from the previous time step and p c denotes 
the coarse-scale nodal values of the discrete CG solution. We also note that the operator matrix may be 
analogously used in order to project coarse scale solutions onto the fine grid 

pu+i = ( RP f p n+l 


4.2 Displacement Solve 

We now suppose that we have solved for the fine-grid pressure p n+1 by the GMsFEM pressure solve in the 


previous section. We must now solve the mechanics equations (7a). Since the construction of the multiscale 


offline space remains very similar in this setting, we will be a bit more brief on its construction. Recall, for 
discretization of the displacements equation we rewrite equation as follows 


Au 


n+1 


= F U 


(15) 


where F u = —aGp n+l . The corresponding continuous Galerkin (CG) formulation for displacements equa¬ 
tions is given by: 

a{u n+1 ,v) = (f u ,v), for all v £ V 0 s, (16) 
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where u(x,t) = Si fc K WK'( x )> where are fine-scale basis functions, and we construct the multiscale 
offline space V 0 g. 

For construction of multiscale basis functions for displacements we use similar algorithm that we used 
for pressure. For construction of a snapshot space F s ^ ap we solve following problem in oj 

a(v>j P ,v)=0 mu, 

ip“' snap = Sj( X ), onduj. 

Let Li be the number of functions in the snapshot space in the region oj, and define 

Knap = span{<^ nap : 1 <j<k}, 

for each coarse subdomain oj. Note we are using the same notation but with different harmonic extensions. 
We denote the corresponding matrix of snapshot functions, again with similar notation, to be 

Knap =[< aP ,...,< aP ]. 

Again, we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decom¬ 
position. We solve the eigenvalue problem in the space of snapshots 

A off $f = A (18) 

where where 

^° ff = (Knapf^Knap, N° S = (Knapf^Knap, 

where A and N denote fine scale matrices 


a = 

1 


J ^ ( y^(<Pm ) : £ (V>J + Adiv (^m) ' div Oj)> 


and 


N-mn — / ( A + • lf n . 

J n 


Here, ip, are fine-scale basis functions. 

We then choose the smallest N^g U eigenvalues from Eq. (18) and form the corresponding eigenvectors in 
the space of snapshots by setting 

3 =1 


snap 
°3 ’ 


for k = 1,..., N^g U , where are the coordinates of the vector ip'jf. We denote the span of this reduced 
space as V£ s . 

For construciton of multiscale partition of unity functions for the mechanics solve, we proceed as before 
and solve for all I\ £ oj 


a(€i,v) = 0 in K, (19) 

£i = gi on dK. (20) 

Here gi is a continuous function on K and is linear on each edge of dK. Finally, we multiply the partition 
of unity functions by the eigenfunctions in the offline space to construct the resulting basis functions 

<Pi,k = f° r 1 < i < N and 1 < k < N^’ u , (21) 

where N“g U denotes the number of offline eigenvectors that are chosen for each coarse node i. Next, we 
define the spectral multiscale space as 

Vos = spanjy^fc : 1 < i < N and 1 < k < N^}. (22) 



Using a single index notation, we may write V 0 s = spanjt/jj}^, where = J2iLi ^off U denotes the total 
number of basis functions in the space V£g, for all i = 1,..., N. 

And after construction V 0 ff we denote the matrix 

R U = [+ 1 , • ■ ■ ) t PN?] T , 


where ipi are used to denote the nodal values of each basis function defined on the fine grid. Then, the 


variational form in (16) yields the following linear algebraic system 


A c u, 


n+1 


= F r 


(23) 


where A c = R U A(R U ) T , F c = R U F U and 


u n+i = i 


5 Numerical Examples 


In this section, we present numerical examples to demonstrate the performance of the GMsFEM for comput¬ 
ing the solution of the poroelasticity problem in heterogenous domains. Although we presented the algorithm 
in the fixed-stress splitting, we are able to apply the same offline spaces (Q 0 s, I4ff) as their construction re¬ 


mains the same in the fully coupled setting. However, in the coupled setting the equations (141 and (23) will 
no longer be decoupled and must be solved simultaneously. 

We will implement a single complicated geometry with contrasting parameter values. We provide two 
cases one with lower contrast in elastic properties and another with higher contrast. We present the algorithm 
applied to these heterogenous coefficients in both the fully coupled and fixed stress time splittings. We give 
the errors with varying multiscale basis functions and over time. We then will apply the GMsFEM method 
with oversampling and with snapshots with randomized boundary conditions to obtain good accuracy, while 
having to solve fewer snapshot solutions. The effects of higher contrast in properties will also be discussed. 


5.1 GMsFEM Implementation 

First, we take the computational domain H as a unit square [0,1] 2 , and set the source term / = 0 in ([l]). 
We utilize heterogeneous coefficients that have different values in two subdomains. We denote each region 
as subdomain 1 and 2, and use following coefficients: for the Biot modulus we take M-\ = 1.0, M 2 = 10 
and for permeability ki = 10~ 3 , k 2 = 1 in the two regions. For fluid viscosity we take v = 1 and fluid-solid 
coupling constant a = 0.9. For the elastic properties, we present results for two test cases. In Case 1, 
the elastic modulus is given by E\ = 10, E 2 = 1 in each respective subdomain and in Case 2, we have 
E\ = 10, E 2 = 10 -3 . The Poisson’s ratio is 77 = 0.22, and these can be related to the parameters //,; and A 
for i = 1,2, via the relation 

- Ei \ - EiT] 

2(1+7?)’ * (1 + f?)(l- 277)’ 

in each subdomain. The subdomains for coefficients shown in Fig. [2j where the background media in red is 
the subdomain 1, and isolated particles and strips in blue are the subdomain 2. 

As we have chosen / = 0 we must use boundary conditions to force flow and mechanics. In these tests, 
we use following boundary conditions: 


and 


P = Pi, xGT t , P = Po, ieffl, 



x €T l u r#, 


u x = 0, 


du y 

dy 


= 0 , 


x g r L , 


du x 

dx 


0, u y = 0, x G Tg, 
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Figure 2: Coefficients subdomains. Red is the subdomain 1 and blue is the subdomain 2 


and finally, 


du x 

dx 


= 0 , 


8Uy 

dy 


= 0 , 


x £ r T u r R . 


Here T R and are left and right boundaries, IV and T b are top and bottom boundaries respectively. We 
set po = 0 and pi = 1 to drive the flow, and thus, the mechanics. 




Figure 3: Coarse and fine grids 


In Fig. [3]we show the coarse and fine grids. The coarse grid consists of 36 nodes and 50 triangle cells, and 
the fine mesh consists of 3721 nodes and 7200 triangle cells. The number of time steps is 20 and the maximal 
time being set at T max = 100. As an initial condition for pressure we use p = Po- The reference solution 
computed by using a standard FEM (linear basis functions for pressure and displacements) on the fine grid 
and using a fully coupled scheme. The pressure and the displacement fields for Case 1 on the fine-grid are 
presented on the left column of Fig. 0-H 

We test the fully coupled and fixed-stress splitting schemes. The errors will be measured in weighted L 2 
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and weighted H 1 norm and semi-norm for pressure 


E pll_L 2 (n) 


-p\h 1 {Q) 


J n ^ p f ~ Pmsfdxj 

J grad (p f - p ms ), grad(p/ - p ms )^ dx'j , 


and for displacements 


^Hl 2 (0) — "h (^/ ^ms)dx^ 

/ (a(u f - u ms ),e(uf - u ms ))dx 
Jn 


1/2 


l e “l//i(n) — 


Here (uf,pf) and ( u m s,Pms ) are fine-scale and coarse-scale using GMsFEM solutions, respectively for pres¬ 
sure and displacements. 

Recall, we will use a few multiscale basis functions per each coarse node w.;, and these number of coarse 
basis defines the problem size (dimension of offline spaces, Q 0 g and V 0 g). We suppose that in each patch w, 
we take the same number of multiscale basis functions for pressure, N p s = N^’ p , for i = 1, ■ ■ ■ ,1V. Similarly 
for displacements we take lV“ ff = N^' u , for i = 1, • ■ • , N. Varying the basis functions in both pressure and 
displacement multiscale spaces we recorded the errors at the final times. 

In Tables [l] and [ 2 J we present the weighted L 2 and H 1 errors for Case 1 and Case 2 of the coefficients in 
geometry Fig. fusing the fully coupled scheme. We compare these to a fine-scale solution space with dimen¬ 
sion 11163. In these tables, N p s and IV“ ff are number of multiscale basis functions for each neighborhoods, 
the second column show the dimension of the offline space, the next two columns present the weighted L 2 
and H 1 errors for pressure and last two columns show the weighted L 2 and H 1 errors for displacements. We 
see that the errors in pressure remain similar in both cases because the permeability parameters remain the 
same and the change is in elastic properties between scenarios. In Case 2, pictured in Table [2j we see great 
errors in displacements throughout when compared to Case 1 in Table [T] because the elastic properties in 
Case 2 have several orders of higher contrast. 

In a similar setting, we consider the fixed-stress splitting. For Case 1 we present the results in Table [3j 
the errors are very similar compared to the corresponding fully coupled scheme. This may be because we are 
comparing a fine-scale fully coupled scheme to a multiscale fully coupled scheme and similarly, a fine-scale 
splitting scheme to a multiscale splitting scheme and the errors do not differ very much between the two 
schemes here. For Case 2 we present the errors in Table [4] and again see that the errors are higher when 
compared to the lower contrast scenario. Comparing these results with the Case 2 using the fully coupled 
scheme, presented in Table [2| we see that both the pressure errors and displacement errors are much greater 
in this sequential coupling. This disparity is particularly striking when few multiscale basis functions are 
used. 

We also include plots over time of the error with respect to number of basis functions used. We present 
the results from the fully coupled scheme. In Fig. [6] and [ 7 ] we show errors over time for N 0 g = N p s = 7V“ ff = 
4, 8,12, and 16 multiscale basis functions for each Ui Thus, the dimensions of offline spaces are 432, 864, 1296 
and 1728, respectively. We observe that errors decrease as we increase the dimension of the offline space as 
expected. We observe the errors in Fig. [6]are generally better than the errors Fig. [TJ again, due to the lower 
contrast in Case 1. We see that in both cases most of the error vanished after the use of just 8 multiscale 
basis functions. In general, the error remains stable in time with a slight decrease over time. 


5.2 GMsFEM with Randomized Oversampling 

In this section we consider the oversampling randomized algorithm proposed in [S] . In this algorithm, instead 
of solving harmonic extensions ([9] and © for each fine grid node on the boundary, we solve a small number 
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*5f 

dim(Q off , V oS ) 

Pressure errors 
L 2 H 1 

Displacements errors 
L 2 H 1 



N u 

iV off 

= 2 



2 

216 

0.06 

0.08 

0.06 

0.13 



N u 

iV off 

= 4 



2 

360 

0.06 

0.08 

0.06 

0.12 

4 

432 

0.01 

0.01 

0.04 

0.11 



N u 

iV off 

= 8 



2 

648 

0.06 

0.08 

0.02 

0.06 

4 

720 

0.01 

0.01 

0.01 

0.03 

8 

864 

0.0003 

0.002 

0.002 

0.03 



N u - 
iV off ■ 

= 12 



2 

936 

0.06 

0.08 

0.02 

0.05 

4 

1008 

0.01 

0.01 

0.01 

0.02 

8 

1152 

0.0003 

0.002 

0.0009 

0.01 

12 

1296 

0.0001 

0.001 

0.0009 

0.01 



N u - 
iV off ■ 

= 16 



2 

1224 

0.06 

0.08 

0.02 

0.05 

4 

1296 

0.01 

0.01 

0.01 

0.01 

8 

1440 

0.0003 

0.002 

0.0008 

0.01 

12 

1584 

0.0001 

0.001 

0.0007 

0.01 

16 

1728 

0.0001 

0.0007 

0.0007 

0.01 


Table 1: Numerical results for Case 1 using the fully coupled scheme. 



dim(Q 0 ff, V 0 g) 

Pressure errors 
L 2 H 1 

Displacements errors 
L 2 H 1 



N u 

= 2 



2 

216 

0.06 

0.08 

0.25 

0.26 



N u 

iV off 

= 4 



2 

360 

0.06 

0.08 

0.22 

0.24 

4 

432 

0.02 

0.01 

0.19 

0.24 



N U 

iV off 

= 8 



2 

648 

0.06 

0.08 

0.08 

0.13 

4 

720 

0.02 

0.01 

0.01 

0.08 

8 

864 

0.001 

0.002 

0.01 

0.08 



Afu . 
iV off ' 

= 12 



2 

936 

0.06 

0.08 

0.07 

0.11 

4 

1008 

0.02 

0.01 

0.02 

0.04 

8 

1152 

0.0003 

0.002 

0.004 

0.03 

12 

1296 

0.0001 

0.001 

0.004 

0.03 



N u - 

= 16 



2 

1224 

0.06 

0.08 

0.07 

0.11 

4 

1296 

0.02 

0.01 

0.02 

0.03 

8 

1440 

0.0003 

0.002 

0.001 

0.02 

12 

1584 

0.0001 

0.001 

0.001 

0.02 

16 

1728 

0.0001 

0.0006 

0.001 

0.02 


Table 2: Numerical results for Case 2 using the fully coupled scheme. 
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N P oS 

dim(Q off , V 0 g) 

Pressure errors 
L 2 H 1 

Displacements errors 
L 2 H 1 



N u 

iV off 

= 2 



2 

216 

0.06 

0.08 

0.06 

0.13 



N u 

iV off 

= 4 



2 

360 

0.06 

0.08 

0.06 

0.12 

4 

432 

0.01 

0.01 

0.04 

0.11 



N u 

iV off 

= 8 



2 

648 

0.06 

0.08 

0.02 

0.06 

4 

720 

0.01 

0.01 

0.01 

0.03 

8 

864 

0.0003 

0.002 

0.002 

0.03 



N u - 
iV off ■ 

= 12 



2 

936 

0.06 

0.08 

0.02 

0.05 

4 

1008 

0.01 

0.01 

0.01 

0.02 

8 

1152 

0.0003 

0.002 

0.0009 

0.01 

12 

1296 

0.0001 

0.001 

0.0009 

0.01 



N u - 

= 16 



2 

1224 

0.06 

0.08 

0.02 

0.05 

4 

1296 

0.01 

0.01 

0.01 

0.01 

8 

1440 

0.0003 

0.002 

0.0008 

0.01 

12 

1584 

0.0001 

0.001 

0.0007 

0.01 

16 

1728 

0.0001 

0.0007 

0.0007 

0.01 

’able 3: 

Numerical results for Case 1 using the fixed-stress scheme 


dim(Q off , V 0 s) 

Pressure errors 
L 2 H 1 

Displacements errors 

L 2 H 1 



N u 

iV off 

= 2 



2 

216 

0.30 

0.26 

0.45 

0.46 



N u 

- /V off 

= 4 



2 

360 

0.30 

0.26 

0.42 

0.45 

4 

432 

0.01 

0.01 

0.33 

0.38 



N u 

iV off 

= 8 



2 

648 

0.30 

0.25 

0.36 

0.48 

4 

720 

0.006 

0.01 

0.04 

0.15 

8 

864 

0.001 

0.006 

0.04 

0.15 



Afu - 
iV off ■ 

= 12 



2 

936 

0.30 

0.25 

0.37 

0.50 

4 

1008 

0.006 

0.01 

0.007 

0.06 

8 

1152 

0.002 

0.006 

0.007 

0.06 

12 

1296 

0.001 

0.004 

0.007 

0.06 



N u - 

= 16 



2 

1224 

0.30 

0.25 

0.38 

0.50 

4 

1296 

0.006 

0.01 

0.003 

0.03 

8 

1440 

0.002 

0.006 

0.002 

0.02 

12 

1584 

0.001 

0.004 

0.002 

0.02 

16 

1728 

0.0009 

0.003 

0.002 

0.02 


Table 4: Numerical results for Case 2 using the fixed-stress scheme. 
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of harmonic extension local problems with random boundary conditions. More precisely, we let 

/Ui.snap _ o 

tpj = r j> x € duii, 

where rj are independent identical distributed standard Gaussian random vectors on the fine grid nodes of 
the boundary. The advantage of this algorithm lies in the fact that a much fewer number of snapshot basis 
functions are calculated, while maintaining accuracy. In addition, we will use an oversampling strategy. This 
is done to reduce the mismatching effects of boundary conditions imposed artificially in the construction 
of snapshot basis functions. We will denote the extended coarse grid neighborhood for t = 1,2,..., by 
w+ = 0 Ji +1. Here for example, + 1, would mean the coarse grid neighborhood plus all 1 layer of 

adjacent fine grida of Wj, and so on. 

In Fig. [8] and Fig. [f)J we show the weighted L 2 and H 1 errors over time for Case 1 and 2 using 
the randomized GMsFEM with oversampling using different numbers of multiscale basis functions. The 
oversampled region + 4 is chosen, that is, the oversampled region contains an extra 4 fine grid cells 

layers around Wi. Here, we use only the fully coupled scheme. We use a snapshot ratio of 36% between 
the standard number of snapshots and the randomized algorithm. Comparing results from Fig. [8j the 
randomized algorithm, to Fig. [6j the standard GMsFEM, we observe that the randomized algorithm is 
slightly less accurate but at the advantage of having less snapshot solutions required. 

In Table [ 5 ] and [g] we investigate the effect of the oversampling uif = Ui + t as we increase the number 
of fine grid extensions for t = 0,2,4 and 6. We present the data of the randomized snapshots for last time 
step. We see that oversampling can help to improve the results initially, but the improvements level off as 
large oversampling domains do not give significant improvement in the solution accuracy. Again the effects 
of the high contrast of Case 2 can be seen in the data as the oversampling performs slightly worse than in 
the lower contrast regime. 


N 0 g 

pressure errors 
L 2 H 1 

displacements errors 
L 2 H 1 


without oversampling, 

W+ = Uli 

4 

0.05 

0.04 

0.16 

0.22 

8 

0.05 

0.03 

0.16 

0.21 

12 

0.02 

0.02 

0.16 

0.21 

16 

0.004 

0.01 

0.15 

0.21 


with oversampling, w z + 

= 2 

4 

0.05 

0.03 

0.14 

0.21 

8 

0.04 

0.03 

0.12 

0.19 

12 

0.007 

0.01 

0.09 

0.17 

16 

0.002 

0.009 

0.08 

0.16 


with oversampling, uj^ 

= Wi + 4 

4 

0.05 

0.03 

0.09 

0.17 

8 

0.04 

0.02 

0.06 

0.14 

12 

0.006 

0.01 

0.04 

0.11 

16 

0.001 

0.008 

0.02 

0.08 


with oversampling, w z + 

= (jOi + 6 

4 

0.05 

0.03 

0.09 

0.17 

8 

0.04 

0.02 

0.06 

0.13 

12 

0.009 

0.01 

0.02 

0.09 

16 

0.002 

0.007 

0.02 

0.07 


Table 5: Numerical tests for Case 1 using randomized GMsFEM with and without oversampling for N a g = 
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pressure errors 

displacements errors 1 

N 

L 2 

H 1 

L 2 

H 1 


without oversampling, 

= UJi 

4 

0.05 

0.04 

0.36 

0.31 

8 

0.05 

0.03 

0.35 

0.31 

12 

0.02 

0.02 

0.34 

0.31 

16 

0.006 

0.01 

0.34 

0.31 


with oversampling, uj^ 

= uji 2 

4 

0.05 

0.03 

0.33 

0.31 

8 

0.04 

0.03 

0.30 

0.30 

12 

0.01 

0.01 

0.25 

0.27 

16 

0.009 

0.009 

0.22 

0.25 


with oversampling, 

= LOi + 4 

4 

0.05 

0.03 

0.28 

0.29 

8 

0.03 

0.02 

0.19 

0.24 

12 

0.007 

0.01 

0.11 

0.20 

16 

0.002 

0.008 

0.07 

0.15 


with oversampling, 

= Cc^ + 6 

4 

0.05 

0.03 

0.24 

0.27 

8 

0.04 

0.02 

0.14 

0.22 

12 

0.01 

0.01 

0.07 

0.17 

16 

0.002 

0.007 

0.06 

0.14 


Table 6: Numerical tests for Case 2 using randomized GMsFEM with and without oversampling for N () ff 
^off = Ks 
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6 Conclusion 


Simulating poroelasticity is difficult due the complex heterogeneities and because of the complexity of grid- 
ding the flow and mechanics regimes in such media. Therefore, in this paper we developed a Generalized 
Multiscale Finite Element Method for a linear poroelastic media. We first presented the general poroe¬ 
lasticity framework of Biot and its subsequent solution by fixed stress time splitting methods. Although 
fully coupled schemes are considered numerically, this splitting lays the framework for the application of the 
GMsFEM to the decoupled poroelastic equations. We then outline the construction of the multiscale spaces 
in both fluid and mechanics regimes. The algorithm is then implemented on a single geometry with two 
different cases of elastic parameters. We show the errors relative to the fine scale solution over time and with 
varying multiscale basis functions. Finally, we implemented oversampling strategies and randomized bound¬ 
ary conditions when solving for the snapshot space. As in cases of reservoir compaction, the permeability 
may depend on pressure resulting in a nonlinear relation. In future studies, we will develop a GMsFEM for 
such nonlinear poroelastic problems. 
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Fine 


Coarse 



Figure 4: The fine-scale and coarse-scale solutions of the pressure distribution for T = 10 and 100 (from 
top to bottom) for case 1. The dimension of the fine-scale solution is 11163 and the dimension of the coarse 
space is 864. 
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Fine Coarse 



Figure 5: The fine-scale and coarse-scale solutions of the displacements u x and u y for case 1. The dimension 
of the fine-scale solution is 11163 and the dimension of the coarse space is 864. 
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time 


Figure 6: Weighted L 2 are on the top and H 1 are on the bottom. Errors for pressure are on the left and 
displacements are on the right for Case 1 using the fully coupled scheme. 
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Figure 7: Weighted L 2 are on the top and H 1 are on the bottom. Errors for pressure are on the left and 
displacements are on the right for Case 2 using the fully coupled scheme. 
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Figure 8: Weighted L 2 are on the top and H l are on the bottom. Errors for pressure are on the left and 
displacements are on the right for Case 1 using randomized GMsFEM with oversampling, uf = + 4. 
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Figure 9: Weighted L 2 are on the top and H l are on the bottom. Errors for pressure are on the left and 
displacements are on the right for Case 2 using randomized GMsFEM with oversampling, uf = + 4. 
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